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Abstract We present a method for determining surface flows from solar im¬ 
ages based upon optical flow techniques. We apply the method to sets of images 
obtained by a variety of solar imagers to assess its performance. The opf lowSd 
procedure is shown to extract accurate velocity estimates when provided per¬ 
fect test data and quickly generates results consistent with completely distinct 
methods when applied on global scales. We also validate it in detail by com¬ 
paring it to an established method when applied to high-resolution datasets 
and find that it provides comparable results without the need to tune, filter 
or otherwise preprocess the images before its application. 


1 Introduction 

The most powerful events observed in the solar system are the result of convec¬ 
tive dynamics in the outer layers of the Sun. Researchers trying to understand 
these dynamics need tools to measure and study these turbulent motions and 
their interaction with the local magnetic field. Several methods for deducing 
the flows in the solar photosphere (the visible surface of the sun) have been 
developed over the years. Many of these have been based upon some form of 
local correlation tracking (LCT) pillTS] where pairs of successive images of 
the photosphere are broken into sub-images which are then shifted relative to 
each other to find an optimal relative shift. This shift is then associated with 
the local mean velocity of the flow. These methods appear to work well on 
high-resolution images collected from both ground-based [2] and space-based 
mm observatories. There have been a few attempts to assess and compare 
the performance of such methods over the past twenty years including those of 
Hurlburt et al. [5], Welsch et aZ.[2S], Chae and Sakurai [J] and most recently, 
Verma, Steffen and Denker |211. These have used a variety of data input for 
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their assessments, from those derived from simulations of MHD flows in the 
photosphere, to creating image sets from solar images by applying known dis¬ 
tortions, to direct comparisons with real solar data. All tested methods gave 
consistent results. However the best method was somewhat dependent on the 
test and the choice of the respective parameters for each method. 

Alternative methods have been developed in other fields with similar goals. 
Machine vision researchers developed optical flow methods for deducing the 
relative motions of objects in digital images wm and atmospheric scientists 
developed various methods for deducing winds on Earth from cloud motions 
in satellite images mm- Experimentalists in other branches of fluid dynam¬ 
ics, including blood flow measurements deduced from x-ray images m and 
flows with various tracer particles in suspensions [27] . have explored similar 
methods. One difference between the typical machine vision problem and that 
of deducing fluid velocities is that the former seek motions of discrete ob¬ 
jects while the latter seeks motions of continuous flows. Applying optical flow 
methods such as those described by Jahne m on high-resolution solar images 
typically underestimate the velocity of the flows. In part this is due to the 
simple spatial averaging used in deriving the velocity. 

Hurlburt U presented a method which does not make use of such spatial 
averaging. Instead the flows are assumed to be smooth and continuous - being 
represented by a truncated Fourier series. Here we present a detailed descrip¬ 
tion of the method and assess it using previously developed tests and compare 
it to local correlation tracking methods. 


2 Method 

The basic assumption is that the visible pattern observed in the fluid, as 
measured by the local intensity I, will be advected by the velocity field v and 
hence should satisfy the equation 

§+v.V7 = 0 (1) 

In the case of solar physics the pattern is typically formed by convective mo¬ 
tions in the photosphere, which are clearly visible in white-light images and 
which appear to be advected by larger scale flow fields. Images are collected 
frequently relative to the flow speeds, such that the displacement of the pat¬ 
tern between any sequential pair of images is less than a pixel. The problem 
is to determine v from a time sequence of two dimensional images I{x,y,t) 
in the presence of measurement noise and other “noise” sources, such as the 
acoustic oscillations present in the solar atmosphere or missing frames due to 
data dropouts. Using equation Q we can seek the best fit velocity field v/ 
using least squares. First we form the merit function of the fit for the full 
dataset /(x, y, t). 
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which we seek to minimize. Here the sums are taken over the discrete pixel 
and frame coordinates for x,y and t. If the velocity field v is a continuous 
field, we can express the fit velocity Vf as a Fourier series 

V/ = (3) 

where x and y are unit vectors, and /3y are complex amplitudes and and 
Ny are the number of Fourier modes retained in the expansion. Substituting 
this into equation ([^ and differentiating with respect to au and /3ki results 
in the system of equations 

StSrEy{^+Vf ■ ^ Q 

This can be reorganized to form the complex system of equations 

= Rh, = Rl^ ( 6 ) 

where 

nx _ _y y y f 91 dl \ 2-xl(kx! X+ly j Y) 
nV _ _y y y f -2-iTl{kx/X+ly/Y) 

Mfyh = E^SyEt(^^^ ^-2^H{kYi)xlX+{lYi)ylY) 
yrxy _ y y y f -2-iTl{{k+i)x/X + {l+j)y/Y) 

= E.EyEti^^'^ g-2.i((/c+.wx+(/+jWr) 

These matrices consist of discrete Fourier transforms of the various, time- 
averaged products of the spatial and temporal derivatives of the image. The 
matrices of the complex linear system ([^ for the spectral amplitudes atj , Pij 
can be combined to form a single hermitian matrix of size {SN^ x Ny)'^. 

This method has been implemented in an IDL routine opf low3d and 
is available as part of the SolarSoft environment [S]. The time derivative is 
evaluated using finite differencing between sequential images while the spatial 
derivatives is evaluated using 4th order finite differences on the average of the 
two images used for the time derivative. The matrices are then computed for 
the entire time-space cube /(x, i/, t) and the system is solved. 

Solving the system using direct methods can quickly become expensive, 
scaling as {N^ x Ny)^. The method requires 70 seconds on an 2013-vintage 
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Fig. 1 Comparison between a derived velocity field and input field used to distort a solar 
image (background). The input (derived) velocity field is displayed with black (white) arrows 
whose areas are proportional to the magnitude of the local velocity. The relative scale of the 
white arrows has been reduced slightly for aid in comparison. 


iMac workstation (with 3.4GHz processor and 32GB of memory) to obtain 
the solutions for N^. = Vj, = 24 on a 1024 x 1024 image. This is partially offset 
by the fact that the method requires no preprocessing or filtering and that 
it fits flows over many instances in time in one go. The matrices in Q can 
also be reused in subsequent calculations with minimal additional expense. 
Performance could be further optimized by taking advantage of the matrix 
structure in equation 0 which has a blocked-Toeplitz-Toeplitz-block (BTTB) 
structure [3]. 


3 Evaluation and comparisons 

As a first test of the method, we take the simulated observations developed by 
Hurlburt et al. [6]. Using a sixth-order accurate numerical scheme [9] they took 
a single intensity image Iq of solar granulation and evolved it with equation 0 
with a steady velocity field. They then degraded and resampled the image to 
represent the expected resolution of the Michelson Doppler Imager (MDI) on 
the Solar and Heliospheric Observatory m- Since there is no source of noise 
and the imposed flows are themselves based on Fourier modes, we expect and 
observe that the opf low3d method can recover the flow with a high accuracy. 
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Latitude 

Fig. 2 The differential rotation of the sun as determined by applying the spectral optical 
flow technique on one hour of full-disk MDI continuum images (solid) and published best 
fit from Doppler measurements (dotted). The former was calculated with Nx =4, Ny = 28 
using sixty 256 x 896 pixel images. The two agree within the noise level of the supergranular 
flow field. 


The results for a case where N^: = Ny = 4 on a 140 x 140 pixel image is 
displayed in figure along with the known input velocty field and a sample 
image. The two sets of arrows, corresponding to the known (black) and derived 
(white) velocity fields are almost perfectly correlated, both in direction and 
magnitude. 


With this basic validation of the method on perfect data, we turn to com¬ 
parisons with other methods and operate on real solar data and then consider 
a detailed error analysis. First we compare the results of applying this method 
to other large-scale, full-disk measurements. Figure [^displays the zonal (E-W) 
component of the solar velocity along the central meridian of the Sun as a func¬ 
tion of solar latitude derived from one hour of MDI m data. We include the 
corresponding measure based upon fits to Doppler measurements m- Aside 
from the departures induced by sampling errors of the supergranular flows in 
this short time, the two curves agree very well. 
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Fig. 3 The horizontal velocity profile through the middle of the field of view. The magnitude 
of the simulated flow |V| is displayed in black and the fitted velocity |vy |, and error, |v —vy | 
in red and green respectively. It is clear that error drops rapidly away from the boundaries, 
within half a wavelength of the truncated mode. 


4 Error Assessment 

There are several factors that may contribute to errors in the velocity estimate 
provided by opflowSd. These can be broken into three classes: systematic errors 
introduced by the choice of velocity representation, errors due to image quality 
and errors introduced by physical effects in the solar atmosphere. As a first 
step, we take what is currently the most consistent and stable images available, 
using data obtained by the Helioseismic and Magnetic Imager (HMI) on SDO 
m- We then subject them to a variety of controlled tests to address the 
first two classes of error. This is the best-case scenario for studies of the solar 
photosphere: a stable imager observing ’’quiet” sun (where magnetic effects 
are negligible). The following section considers a more complex situation of 
observing magnetic regions with a less-stable imager. 

One hour of Level-1 HMI continuum images consisting of 80 individual 
frames were used for this study beginning at 2010-10-26T08:29:00. A set of 
1024x1024 pixel sub-images were extracted centered on a coronal hole iden¬ 
tified in the Heliophysics Event Knowledgebase [10] (herein referred to as 
C2010). 

The use of a truncated Fourier representation for the velocity field is a 
common practice in fluid dynamic investigations. However it has two well- 
known issues that must be considered: it imposes a periodic structure to the 
flows and may introduce aliasing or other errors due to truncation. To assess 
these effects, we use the first frame from the C2010 and generate ten artificial 
images from it using equation Q with a known, hexagonal flow field, in a 
similar approach to that of Hurlburt . 

Our implementation uses a fast Fourier transform and can be sensitive to 
discontinuities at frame boundaries. While the Fourier representation forces 
the velocity to match across opposing boundaries, the application of a least 
squares fit works to confine such effects to near the edges. The resulting spikes 
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in the residual |v — v^| decay rapidly away from frame boundaries (Figure]^. 
Thus, avoiding pixels near the boundaries is the practical means for avoiding 
these errors. 

Aspects of image quality that may influence our fit include large displace¬ 
ments arising from insufficient sampling rates (which may introduce ambiguity 
into the possible solution), image noise and missing data. To assess the impact 
of image quality, we adapt the approach from the previous section. The first 
frame of the C2010 is advected to generate ten frames with a known velocity 
field V composed of 22 by 22 randomized Fourier modes. A set of these data 
cubes is generated with a range of magnitudes for V. We find the accuracy of 
the fit starts to break down when the RMS velocity exceeds 0.4 pixels/frame. 
This result applies to the images with repeating, high-frequency patterns such 
as solar granulation seen in C2010. In contrast, larger-scale patterns could 
allow larger unambiguous motions between frames. 

The effect of truncation can be assessed by evaluating this same test case 
with an increasingly large number of modes (Nx,Ny) in equation (|^, from 
Nx,Ny = 4 to 16. We find that the overall trend is constant, but with higher 
values being more sensitive to the effects of noise as the effective sample size 
of the fit decreases. Thus the selection of the number of modes should take 
this trade off into account. A rule of thumb would be that < Nf, where 
Nf is the number of features (e.g. granules) required to span the image. 

To assess the impact of noise on velocity estimate, we first measure the in¬ 
herent noise in the synthetic datacube used above by setting v/ = 0 in equation 
([^ to give x(0) = RMS(^) = 700 counts per frame. x(0) is used here to rep¬ 
resent the original variation that is reduced by fitting V, and because it scales 
with image contrast. Next, we generate a sequence of datacubes by adding 
increasing levels of Gaussian white noise. The datacube with added noise stan¬ 
dard deviation a has velocity estimate v/(cr). By comparing RMS(v/((t) — V) 
with RMS(V), we can gain some insight into the sensitivity of the method 
to noise. We find that for a maximum relative error of 1%, the maximum 
additional noise must have a < 200. Comparing this value of a to the in¬ 
herent signal noise x(0)j demonstrates that high noise images can still yield 
reliable velocity estimates. This robustness against measurement noise most 
likely results from averaging over many pixels. 

Finally there are features in the solar atmosphere that may impact the 
performance of any method of velocity estimation. These include the presence 
of strong acoustic modes (known as five-minute oscillations) which generate 
a relatively-smooth, but random intensity fluctuation in solar images; Limb- 
darkening, which introduces fixed, large-scale intensity gradients due to line- 
of-sight effects; and strong magnetic features that may distort the intensity 
patterns in non-obvious ways. In exploring these cases, we cannot compare 
our results to a known solution. Instead we must make statistical inferences. 

If we seek flows that persist on time scales significantly larger than five 
minutes, we can assess the effect of acoustic oscillations. We take v/(W j) 
to be the velocity fit for N frames starting at frame j. We examine the con¬ 
vergence of Vf{N,j), as N increases. Using C2010, we divide the 80 frames 
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into sets 1 and 41, consisting of the first and last sets of 40 frames. Since 
the images in this case may possess an overall motion akin to camera motion 
(say Vc(V,j)), we first subtract such motions from v^(40,1) to produce the 
the Euclidean metric RMS(|v/(40,1) — Vc(40,1))|) = 0.18 pixels/frame. Sim¬ 
ilarly, the distance between vy(40,l) and vy(40,41) is the Euclidean metric 
RMS(|v/(40,1) — v/(40,41)|) = 0.019 ppf. This distance is a rough measure 
of precision. Thus, with IV = 40, considerable convergence is apparent. 

To examine the rate of convergence, we can next estimate vy(2,j) for 
j = 1,39. We then compute the distance RMS(v/(2, j) — v/(40,41)). This 
distance, averaged over all j, is 0.19 ppf. Clearly a two frame estimate is 
poor. Comparison of the 2 frame and 40 frame distances, 0.19 ppf and 0.019 
ppf, imply that the rate of convergence is between 1 /-s/fV and 1 /N, which is 
expected for traveling wave patterns in acoustic oscillations. This is consistent 
with previous studies of solar flows using LCT methods m- 

The contrast in the image varies smoothly across each frame in our sample 
due to the limb-darkening effect of the solar atmosphere. Such gradients can 
cause problems for some methods. However in our case, it only changes the 
weights of the sum of the squares in equation ([^, so that the fit is only 
slightly affected. The same logic shows the method is insensitive to image- 
quality problems such as missing frames. 


5 Comparison to other methods 

As a final test, we provide a detailed comparison of our method to that used 
a recent study by Verma & Decker [5S] (hereinafter referred to as V&D). 
They conducted a thorough investigation of horizontal flow fields observed in 
Hinode G-band images using a local correlation tracking (LCT) approach that 
was used in November and Simon US]. 

Following V&D, we selected an hour-long set of G-band images collected 
by the Hinode Solar Optical Telescope (SOT, [23]) on June 4, 2007 between 
14:27 and 15:27UT. In that study, the authors first applied the standard cal¬ 
ibrations to the images and then further pre-processed them by correcting 
for foreshortening, applying a rigid alignment between the images to remove 
spacecraft jitter and solar rotation, and then employed a subsonic filter to 
remove acoustic oscillations. 

We also calibrate the images using the SolarSoft routine fg_prep with 
its default settings. However, we do not apply the other preprocessing steps. 
Other than foreshortening, those corrections effectively remove noise from the 
velocity signal that we are seeking, be it jitter from the spacecraft, bulk motion 
across the field of view or distracting intensity fluctuations. Since the opf low3d 
method has already been shown to address such noise sources, we rely on it 
alone to do so. In addition, we found seven of the 238 images in the sequence 
were missing: rather than attempt to correct for these, we left those images 
blank and left it to opflowSd to deal with the consequences. 
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Fig. 4 A comparison of results of the method used by (a) V&:D (from their Figure 2c) and 
(b) opf low3d for the same sample of Hinode/SOT data shows detailed agreement. Here we 
display the magnitude of the two velocity fields using approximately the same color map 
and scaling from black/dark blue to yellow/red. The flow velocities exhibit the same pattern 
of outward moat flows around sunspots and inflows around plage. 


Figures]^ and [^display a comparison of applying the two methods to the 
same image set. The only free parameters for opflowSd are the number of 
modes used to fit the velocity field, and whether to use a direct or iterative 
solution method: we select a direct solution with 20 modes in each direction. 
This corresponds to an effective pixel size of about 2Mm when compared to 
the Gaussian FWHM used by V&D. 

We correct for foreshortening after the fact by scaling the components of 
the velocity and display the magnitude of the resulting velocity field (less the 
average velocity over the frame) in Figure Strong moat flows are visible 
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Fig. 5 Histograms of the speeds found by (a) V&D (from their Figure 6) and (b) opf low3d 
for the same one-hour sample of Hinode/SOT data also show detailed agreement. The solid 
curves display the normalized histogram of the velocities for the two methods. (The other 
curves in (a) correspond to evaluations over longer time intervals from 2 to 16 hours.) 


in the lower left, as well as converging flows elsewhere in plage areas. Figure 
displays a normalized histogram of flow speeds, which can be compared to 
figure 6 of V&D. In both cases the peak value across the held of view is around 
0.3 km/s. We hnd the overall rms speed to be 0.46km/s, the median to be 0.42 
km/s and the maximum to be 1.68; as compared to 0.44 km/s, 0.40 km/s and 
1.95 km/s respectively. The fact that opf low3d method retains slightly higher 
rms velocities while reducing the extremes suggests that it might both retain 
a higher resolving power while mitigating the inhuences of outliers. 


6 Discussion 

We have described a method for deriving hows from sets of images obtained 
by a variety of solar imagers. The opf low3d procedure has been shown to ex¬ 
tract accurate velocity estimates when provided perfect test data and quickly 
generates results consistent with completely distinct methods when applied on 
global scales. We have also verihed that it agrees in detail with an established 
method when applied to high-resolution datasets - and without the need to 
tune, hlter or otherwise preprocess the images before its application. It is cur- 
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rently running as part of the HEK system to identify regions of solar eruptions 
[8] from data collected by the Atmospheric Imaging Array on SDO [12]. 

Our method has been found to work well on other types of image data, 
including magnetograms, since the only assumptions made are that the mo¬ 
tions displayed in them are reasonably smooth and persistent. It can also be 
combined with other image processing methods to extract motions of spe¬ 
cific features within the field. For instance the motion of the two polarities 
(North/South) in magnetograms could be tracked by thresholding the images 
prior to using opflowSd. Similarly, particular scales could be extracted by 
using high- or low-pass filters. 

With the basic approach established, there are several avenues for improve¬ 
ment. First, we could replace the model equation 0 with a more elaborate 
one, say one that solves the vertical component of the induction equation to 
extract velocities from sets of magnetograms. Second, we could provide a more 
elaborate fitting function, say one that permits a simple time dependence. Fi¬ 
nally one can seek to optimize the method using more sophisticated tools of 
linear algebra. We will explore some of these options in future work. 
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